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ABSTRACT 

In this work, we investigate the ecHpse timing of the polar binary HU Aquarii that has been ob- 
served for almost two decades. Recently, Qian et al. attributed large (0-C) deviations between 
the eclipse ephemeris and observations to a compact system of two massive jovian compan- 
ions. We improve the Keplerian, kinematic model of the Light Travel Time (LTT) effect and 
re-analyse the whole currently available data set. We add almost 60 new, yet unpublished, 
mostly precision light curves obtained using the time high-resolution photo-polarimeter OP- 
TIMA, as well as photometric observations performed at the MONET/N, PIRATE and TCS 
telescopes. We determine new mid-egress times with a mean uncertainty at the level of 1 sec- 
ond or better. We claim that because the observations that currently exist in the literature are 
non-homogeneous with respect to spectral windows (ultraviolet. X-ray, visual, polarimetric 
mode) and the reported mid-egress measurements errors, they may introduce systematics that 
affect orbital fits. Indeed, we find that the published data, when taken literally, cannot be 
explained by any unique solution. Many qualitatively different and best-fit 2-planet config- 
urations, including self-consistent, Newtonian A^-body solutions may be able to explain the 
data. However, using high resolution, precision OPTIMA light curves, we find that the (0-C) 
deviations are best explained by the presence of a single circumbinary companion orbiting at 
a distance of ~ 4.5 AU with a small eccentricity and having ~ 7 Jupiter-masses. This object 
could be the next circumbinary planet detected from the ground, similar to the announced 
companions around close binaries HW Vir, NN Ser, UZ For, DP Leo or SZ Her, and planets 
of this type around Kepler- 16, Kepler-34 and Kepler-35. 

Key words: extrasolar planets — LTT technique — N-body problem — polar — star: HU Aqr 



1 INTRODUCTION 

Magnetic cataclysmic variables (CVs, polars, a.k.a. AM Her stars) 
are interacting close binary systems. They consist of a main se- 
quence red dwarf secondary filling its Roche lobe, and a strongly 
magnetized white dwarf (WD) primary, with typical magnetic field 
values of 10-80 MG (Schwope et al. 2001). The strong magnetic 
field of the primary interacts with the weaker magnetic field of 



the secondary and locks the two stars together. Hence, the syn- 
chronously rotating WD spins at the same rate as the orbital mean 
motion of the binary. Under the gravitational field of the primary, 
material flows from the donor star initially along the binary orbital 
plane, and finally is accreted quasi-radially onto the magnetic poles 
of the WD. The variable HU Aquarii system (hereafter HU Aqr) be- 
longs to this class of CV binaries hosting a strongly magnetic WD 
accompanied by a red dwarf (spectral type M4V) with an orbital 
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period of about 125 minutes. This system is one of the brightest 
polars in the optical domain with visual magnitudes ranging from 
14.6 to 18 (Warner 1995; HelUer 2001), as well as in the X-ray 
energy range. Therefore, it has also been one of the most studied 
systems so far. 

Accreted matter leaving from the red dwarf is initially not 
affected by the magnetic field of the WD. The matter follows a 
ballistic trajectory up to the moment when the WD magnetic field 
begins to dominate. Because the WD magnetosphere extends be- 
yond the Li radius, the plasma stream cannot orbit freely, and thus 
does not form an accretion disk, unlike in other non-magnetic cat- 
aclysmic variables. The accreted matter follows the magnetic field 
lines and forms an accretion spot at the magnetic poles of the WD. 
In many systems, the WD magnetic field is tilted in a such way that 
one magnetic pole is oriented toward the direction of flowing mat- 
ter. Eclipses observed in highly inclined polars provide information 
about the stream geometry. 

According to the most recent work of Schwope et al. (2011) 
the inclination of the binary is ^ 87° ±0.8°. This special geom- 
etry is important for the planetary hypothesis investigated in this 
work. Assuming that a planetary companion (or companions) have 
formed in the circumbinary disc, the inclination constraint removes 
the mass indeterminacy inherent to the eclipse timing method. 

Recently, the HU Aqr system has received much attention in 
the literature. Schwarz et al. (2009) carried out an analysis of the 
light curves of the system and derived mid-egress times of the po- 
lar. They proposed a planetary companion as one possible expla- 
nation of the detected (0-C) variability. Shortly after this work, 
Qian et al. (2011) presented and discussed 10 new light curves in 
the optical domain. These authors confirmed the deviations of the 
observed mid-egress times from a linear or quadratic ephemeris, 
concluding that the large (O-C) residuals may be explained by the 
Light Travel Time [LTT aka Roemer effect, Irwin (1952)] due to 
two jovian-mass planetary companions in orbits with semi-major 
axes of a few AU and a moderate eccentricity of ^ 0.5 for the outer 
planet. The orbit of the inner planet was fixed to be circular. The ra- 
tio of the orbital periods of these massive putative planets would be 
presumably in a low-order 2c: lb mean motion resonance (MMR). 
The latter points to significant mutual interactions between these 
objects which strongly affects the orbital stability of the system. 
Indeed, shortly after that work was published, Horner et al. (201 1) 
performed dynamical analysis of the putative HU Aqr 2-planet sys- 
tem, exploring the parameter space within 3a uncertainty levels of 
the derived Keplerian elements. They found that none of the best- 
fit configurations presented by Qian et al. (201 1) were dynamically 
stable implying that the planetary hypothesis proposed by these au- 
thors is hard to maintain. After a few months, in a new paper, Wit- 
tenmyer et al. (2012) also re-analysed data in Qian et al. (2011) 
confirming that the 2-planet configuration is mathematically con- 
sistent with the observations, but inferred orbits are catastrophically 
unstable over a ~ 10^- 10"^ year time-scale. Furthermore, in a very 
recent paper, Hinse et al. (2012) improved the Keplerian fit models 
of this system by imposing orbital stability constraints on the ob- 
jective function (Xv)^^^- Although these authors were able to find 
a stable 2-planet configuration consistent with the linear ephemeris 
model, orbital parameters of these planets were relatively distant 
from the formal best-fit solution by more than 3a. Because the re- 
sults of extensive dynamical analysis contradict the 2-planet hy- 
pothesis, an alternative explanation of the (O-C) diagrams needs to 
be considered. 

Long term monitoring of HU Aqr shows large variations of 
the accretion rate that could be correlated with a migration of the 



accretion spot. Taking into account the observed changes of the ac- 
cretion geometry during different accretion states, high and inter- 
mediate ones, Schwope et al. (2001) estimated the possible time- 
shift of eclipses to be on the level of 2 seconds, which is still much 
smaller than the deviations between the theoretical ephemeris and 
observed mid-egress moments. These results suggest that the mi- 
gration of the accretion spot cannot be responsible for the (O-C) 
deviations, and we therefore ruled it out. 

The (O-C) variability of HU Aqr could be also attributed to 
other complex astrophysical phenomena in the binary, like the Ap- 
plegate mechanism and/or magnetic braking discussed by Schwarz 
et al. (2009) and Wittenmyer et al. (2012). The timing signal might 
also be affected by non-Gaussian red-noise, which is a well known 
effect present in the precision photometry of transiting planets and 
timing of millisecond pulsars (e.g., Pont et al. 2006; Coles et al. 
201 1). Hence, it should be stressed that we focus here on the plan- 
etary hypothesis, as one of the possible, simple and somehow at- 
tractive explanations of the (O-C) variability. We try to solve "the 
puzzle" of unstable 2-planet models through a new and independent 
analysis of available data, conducted along three basic directions. 

The first one relies on the re-analysis of published data, be- 
cause we found a few inconsistencies in the Hterature. Surprisingly, 
while in the recent paper, Wittenmyer et al. (2012) take into ac- 
count 82 mid-egress points from Schwarz et al. (2009) and Qian 
et al. (2011), this is not the full data set available in the literature 
at that time. In fact, 72 egress times published by Schwarz et al. 
(2009) extend the data set in Schwope et al. (2001) that included 
31 measurements. Although the early data of Schwope et al. (2001) 
spanning cycles 0-22478 overlap with measurements in Schwarz 
et al. (2009) in the time window covering cycles 1319-60097, they 
may be helpful to constrain the best-fit models. Up to now, the full 
list of published observations consists of 1 13 points, including data 
in Qian et al. (201 1). Yet it is not quite obvious whether Qian et al. 
(2011) included measurements in Schwope et al. (2001) in their 
analysis. Hinse et al. (2012) considered the full data set available 
at that time, but in terms of the linear ephemeris LTT model only. 
In this context, a direct comparison of the results in the published 
papers is difficult. 

The second aspect of our study is a new kinematic model of 
the ephemeris that properly approximates orbits of putative com- 
panions in multi-body systems (to the lowest possible order in the 
masses), as compared to the full A/^-body model. The kinematic 
model used in all cited papers refers back to Keplerian parametri- 
sation by Irwin (1952) /or the "one companion" case. That model, 
though commonly used in the literature (e.g., Lee et al. 2011), 
seems nowadays redundant, as it was introduced to quantify a simi- 
larity between the LTT and the radial velocity curves, in a particular 
reference frame with the origin at the center of the two-body LTT 
orbit (instead of the dynamical bary center). Indeed, the recent, al- 
though short history of modeling precision radial velocities teaches 
us that multiple planetary systems should be modeled either us- 
ing kinematic formulation in a proper coordinate frame (e.g. Lee 
& Peale 2003; Gozdziewski et al. 2003), or using the most general 
and accurate full A/^-body model (Laughlin & Chambers 2001). The 
dynamical stability can be further incorporated as an additional, 
implicit observable to the objective function (e.g., Gozdziewski & 
Maciejewski 2001; Gozdziewski et al. 2008). In this work we are 
focused on the kinematic modeling though the self-consistent N- 
body approach was also used to analyse the HU Aqr mid-egress 
times (see Appendix). Our results indicate that the Newtonian 
model may be required for other systems presumably exhibiting 
the LTT effect, indeed. 
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The third and, actually, critical direction of our work, is a 
careful independent analysis of the significantly extended data set 
including already published egress times, and new high-precision 
timing of the egresses obtained with the ultra-fast photometer OP- 
TIMA (Kanbach et al. 2003, 2008), as well as the MONET/N, 
PIRATE and TCS telescopes. We collected almost 60 new egress 
times with superior accuracy at the sub-second level. Moreover, we 
found that the literature data are non-homogeneous, as they come 
from different instruments with different time resolutions, as well 
as working in different spectral windows (from the visual range, 
through the UV, to the X-ray domain) and non/polarimetric modes. 
Taking into account the above mentioned inhomogeneities factors 
and new data, we present the results from a quasi-global optimiza- 
tion of two basic LTT models, leading us to the conclusion that 
the measured (0-C) data of HU Aqr may be best explained by 
a 1 -planet configuration. Simultaneously, it would resolve the 2- 
companions instability paradox in the simplest way. 

This paper is structured as follows. In Sect. 2 we derive 2- 
planet LTT models on the basis of Jacobi coordinates which de- 
scribes kinematic orbits in multiple systems properly, as well as a 
hybrid optimisation algorithm and numerical setup that makes it 
possible to explore the (Xv)^^^ parameter space in a quasi-global 
manner. We also briefly describe the A/^-body formulation of the 
LTT effect. In Sect. 3, we re-analyse the data set published in the 
literature, following the 2-planet hypothesis by Qian et al. (2011) 
and further investigated by Wittenmyer et al. (2012). Two examples 
of highly degenerate best-fit solutions are found. In Sect. 4, possible 
effects of different spectral windows for the light curves and deter- 
mination of egress times are studied. Furthermore, we describe the 
new data set derived with the OPTIMA and other instruments. In 
Sect. 5, we propose the 1 -planet model that best explains the (0-C) 
variability. We briefly discuss the effect of red-noise in Sect. 6 and 
present a summary of our work in the Conclusions, Sect. 7. The 
Appendix contains extensive supplementary material to Sect. 5, in- 
cluding the results of kinematic and A/^-body modeling of 2-planet 
systems, accompanied by the long-term stability tests. 



2 LTT MODEL FOR A 2-PLANET SYSTEM 

We briefly develop the Keplerian model of the LTT signal in the 
three-body configuration, assuming that a compact binary (like 
HU Aqr) has two planetary companions. More technical details 
and a generalization of that model will be published elsewhere 
(Gozdziewski et al., in preparation). We consider the compact bi- 
nary as a single object having the mass of m*, which is reasonable 
in accordance with the extremely short orbital period (~ 125 min) 
of the polar. A single companion, as well as multiple-planet models 
are particular cases of this problem. The key point is that the Ke- 
plerian (or kinematic) model requires special coordinates in order 
to preserve the sense of Keplerian elements as an approximation 
of the exact N-ho&y initial condition. That can be accomplished by 
expressing the dynamics through particular canonical coordinates 
in which the mutual planetary interactions are possibly small with 
respect to the main, "pure" Keplerian part. The barycentric formu- 
lation (Irwin 1952) in fact ignores the interactions which could be 
adequate for low-mass circumbinary objects, but it might fail when 
they have stellar masses as in the SZ Her system (Lee et al. 2011) 
where companions are as massive as 20% of Mq, and can shift the 
system barycenter significantly. The reason for introducing this im- 
proved model is in fact the same as in the precision radial velocities 
analyses (e.g., Lee & Peale 2003; Gozdziewski et al. 2003). 



2.1 Kinematic parametrization of the LTT effect 

One of the well known frames that provides a proper description 
of kinematic orbits in multiple systems is Jacobi coordinates. Let 
us assume that m*, mi and m2 represent the masses of the com- 
pact binary m* and two planets, respectively. Let us also assume 
that the Cartesian coordinates of these objects with respect to the 
three-body barycentre are R*,Ri ,R2, and their Jacobi coordinates 
are denoted by = R*,ri,r2 (see Fig. 1). Here R* is the posi- 
tion of the centre of mass of the binary (CMB) in the barycentric 
frame, and ri , Y2 are position vectors of the planetary companions 
in the Jacobi frame. In this formalism, the barycentric position of 
the binary is: 

R* = -Kiri-K2r2, (1) 

where the mass factor coefficients Ki > 0, K2 > are given by: 
m\ m2 



Ki 



K2 



(2) 



mi + m* mi + m2 + m* 

The coordinate transformation R ^ r is taken from Malhotra 
(1993): 



= R* 



(3) 



r2 = R2 



m*R* +miRi 



mi +m* 

and the inverse transformation is derived from the integral of the 
barycentre: 

R* = -Kiri-K2r2, 

Ri = (1-Ki)ri-K2r2, (4) 
R2 = (1-K2)r2. 

To the first order in the mass-ratio (^ mi 2/^^*), the true A/^-body 
orbit of body /, (/ = 1 , 2) is described through geometric Keplerian 
elements as follows: 



ri{t) = Fi [cos Ei{t) - a] + Qi^Jl-ejsmEi{t), 
where 

Fi = at (1/ cos CO/ + k/ sin co/) , Q/ = a/ (-1/ sin co/ + k/ cos co/) , 
and geometric elements are defined through: 





+ sin Q^i 




+ cos // COS Cli 


1,= 


+ C0S^2/ 


, k,= 


— COS //sin ^2/ 









sin // 



Here, Ei{t) is the eccentric anomaly derived from the Kepler equa- 
tion 

ni{t - Ti) = Ei{t) - eiSmEi{t), 

where nt = 2n/Pi is the mean motion, in accordance with Kepler 
3rd law, njaf = /i/, where Pi is the orbital period of a given object. 
Two tuples {ai,ei,ii,Q.i,(i)i,Ti), i = 1,2, that consist of the semi- 
major axis, eccentricity, inclination, nodal angle, argument of peri- 
centre, and the time of pericentre passage, respectively, are for the 
geometric Keplerian elements. These are related to the Cartesian 
coordinates in the Jacobi frame through the usual two-body formu- 
lae (see, e.g. MorbidelH 2002), with an appropriate mass parameter 
/Hi (see below). 

From Eq. 1, the component of the CMB with respect to the 
system barycentre is: 



Z*(0 = R* . e^ = -KiZi (t) - K2Z2(0: 



(5) 
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to observer 




Figure 1. The geometry of the system. The binary has a total mass 
and because of its short orbital period it can be considered as a point-like 
object accompanied by planets as point-masses. The origin of the coordinate 
system is fixed at the barycentre of the three-body system. The line-of-sight 
is along the z-axis. See the text for more details. 



where is the unit vector along the z-axis of the reference frame, 
directed toward the observer. The signal contribution due to a given 
companion is: 



Zi{t) =aismii 



sin CO/ (cos Ei (t) — + cos co/ yj\ — ej sin Ei (t) 



(6) 



(for planets / = 1 , 2). The Zi (t) are then combined to obtain the (t) 
component of the CMB position vector w.r.t. the system barycentre. 
The LTT signal is then expressed as: 



mi 



x{t) = — Z* = +-^ 

c c \mi -\-m^ 



m2 



mi +m2 +m. 



-Z2 



where c is the speed of light. Note that we used the planetary ver- 
sion of the three-body system, with one dominant mass (m*), hence 
the gravitational Keplerian parameters are: 



jjl=k (mi +m*) 



1^2 ■ 



m*(mi +m2 +m*) 



mi +m* 

consistent with the expansion of the Hamiltonian perturbation for 
the planetary version of the problem (see, e.g., Malhotra 1993), and 
the quantity k denotes the Gauss constant. 

We introduce the signal semi- amplitude factors, Ki and K2 as: 

l\ mi 

J mi +m* 
1\ m2 
) mi +m2 +m=,c 



^1 

^2 



-a\ sm^i, 

-a2 sin/2- 



(7) 
(8) 



Using Eq. 6, the single-planet signal contributions ^, are then given 
by: 



sin CO/ (cosEi(t) — €[) + cos (S^i\j\— e\ sin£/ (t) 



(9) 



In this equation, the set of free orbital parameters is 
(Ki,Pi,ei,(Oi,Ti), i = 1,2, similar to the common kinematic 
radial velocity model. The orbital period P/ and the time of 
pericenter passage are introduced indirectly through the time 
dependence expressed by Ei{t). 



We would like to note here that the contribution of the planet 
as expressed in Irwin's model has an extra term ^/sinco/ that ap- 
pears due to the particular choice of the coordinate system with the 
origin at the center of the binary orbit around the common center 
of mass of the system. It should also be stressed, that no simple 
superposition of kinematic orbits does account for the mutual grav- 
itational interactions directly, but in our formulation, the Keplerian 
elements are the closest to the osculating A/^-body initial condition 
within the kinematic model. 



2.2 The (O-C) formulation 

From Eq. 9, the fit model of the planetary-induced LTT signal is 

T(^7i:i,Pi,^l,COi,ri,i^2,^2,^2,CO2,r2)=ClW + C2(0- 

Now let us assume that the observational data are given through 
eclipse cycle number / (/ = 0, . . . ,N), the date of the eclipse time- 
mark and its uncertainty a/. Then the /-cycle eclipse ephemeris 
with respect to the reference epoch (/ = 0), at time t = ti may be 
written as follows: 

T'ep (/) = ^0 + ^ftin + X (f/ , 7^1 ,2 , A ,2 , ^ 1 ,2 , COi ,2 , 7^1 ,2 ) + "Phy Sics" , 

where Pbin is the orbital period of the binary. It should not be as- 
sumed as known in advance, hence must be fitted, as well as the ini- 
tial epoch ^0 corresponding to cycle number / = 0, simultaneously 
with other parameters of the model. The term coded as "physics" 
contains non-Keplerian effects, such as the period damping or other 
phenomena that may/should be included in the fit model. Here, 
we introduce two instances of such a model. Following Hilditch 
(2001), the linear ephemeris model, as above. 



(O-C) = Tep (/) - ^0 - /Pbin = T (^/ , TTi ,2 , Pi ,2 , ^ 1 ,2 , COi ,2 , Ti ,2 ) 



(10) 



and the quadratic ephemeris model, the simplest, yet non-trivial 
generalization of the polynomial ephemeris (Hilditch 2001) 

(0-C) = rep(/)-?o-/Pbin-p/^=T(^/,7i:i,2,Pl,2,^l,2,COi,2,ri,2).(ll) 

The quantity P in Eq. 1 1 is a factor that describes the binary pe- 
riod damping (change) due to the mass-transfer, magnetic braking, 
gravitational radiation, and/or influence of a very distant compan- 
ion 

P = 2^t)inA)in- 

Let us note that also P should be fitted simultaneously with other 
free parameters of the model. In the rest of this paper, we use a 
common notation in the extrasolar planets literature, that enumer- 
ates the planets by subsequent letters, i.e., "b" = "1", "c" = "2", 
etc., to avoid any confusion. 

2.3 Newtonian model of the LTT effect 

A derivation of the A^-body model of the LTT is basically very 
simple. It requires the computation of the planetary contribution T 
to the (O-C) signal through the numerical integration of the equa- 
tions of motion, a computation of the star barycentric vector R* 
and its -component, in accord with Eq. 5. This formulation ac- 
counts for the mutual interactions between all bodies in the system. 
A serious computational drawback of this model is a significant 
CPU overhead, nevertheless, as we will show in the Appendix, its 
application for systems with massive companions presumably in- 
volved in low order mean motion resonances can be indispensable. 
To solve the equations of motion efficiently, we used the 0DEX2 
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integrator (Hairer et al. 2009) designed for conservative, second or- 
der ordinary differential equations (ODEs). The imposed variable 
time step accuracy preserved the total energy and the angular mo- 
mentum better than 10~^^. In terms of the Newtonian model, the 
planetary bodies are parametrised through the mass mu semi-major 
axis a/, eccentricity et and three Keplerian angles describing the 
orientation of the orbit, for each companion in the system. We also 
assume that the binary is a point mass with the prescribed total 
mass of the binary. Assuming a coplanar configuration (in the N- 
body model the same inclinations are "absorbed" in the planetary 
masses), we have 5 free orbital parameters for each planet, similar 
to the kinematic model. Here, they are then represented as "usual" 
osculating, astrocentric Keplerian elements at a given initial epoch, 
but other types of the osculating elements may be used as well. 

2.4 The optimization method and numerical setup 

Having the egress times measured with a great precision (at the 
1 second level, or even better), the next step is to determine the 
set of primary parameters of the kinematic model, usually with 
the least squares approach, by constructing the reduced (Xv)^^^~ 
squared function 

(X?)l/2 = (Xv)'/'(^l,2,Pl,2,^l,2,C0i,2,ri,2,^0,nin,P), 

and searching for its minimum in the space of the model parame- 
ters. It is well known, however, that the (%v)^^^ function may pos- 
sess many local minima, particularly if the model is not well con- 
strained, as it might be in our case. To seek a global solution, we 
apply a hybrid algorithm that consists of two steps: a quasi-global 
method, the Genetic Algorithm [GA, Charbonneau (1995)] that is 
relatively slow and inaccurate, but makes it possible to find good 
approximations to the second step, a fast local method. Here, we 
use the Levenberg-Marquardt (L-M) algorithm with analytically 
computed derivatives. The idea of the hybrid code comes from 
our earlier works on modelling radial velocity observations (e.g., 
Gozdziewski & Konacki 2004). We used freely available Fortran 
codes of the Genetic Algorithm (PIKAIA^ , by Phil Charbonneau & 
Barry Knapp) and of the L-M method from the well known MIN- 
PACK^ package. 

Once the primary set of the orbital model parameters are de- 
termined in the form of two five-tuples (^/,P/,e/,C0/, 7/), / = 1,2, 
we may also derive inferred Keplerian elements, such as minimal 
planetary masses and semi-major axes, by solving nonlinear equa- 
tions expressing Ki (Eqs. 7, 8) and the 3rd Kepler law in terms of 
the primary model parameters. The inclination has to be held fixed. 
Hence usually one assumes = 90°. Let us underline that while the 
LTT-model (Eq. 9) formulated in the barycenter frame has the same 
mathematical form as in the Jacobi frame, the orbital, geometrical 
(Keplerian) elements in multiple systems should be related to Jaco- 
bian, canonical coordinates. If in the A/^-body numerical integrations 
and stability studies, initial conditions have to be in the form of os- 
culating elements, one should transform these Jacobian elements 
into the Cartesian coordinates w.r.t. the Jacobi frame (e.g., Mor- 
bidelli 2002), and then, if necessary, to the astrocentric or barycen- 
tric coordinates. In this sense, "barycentric" and "Jacobian" two- 
body elements may closely coincide for small, Jovian-mass plan- 
ets. But for more massive companions when the LTT signal is easier 
to detect, or for very compact (resonant) systems, the semi-major 

^ http://www.hao.ucar.edu/modeling/pikaia/pikaia.php 
^ http://www.netlib.org/minpack/ 



axes, masses, Keplerian angles and inferred N-hody initial condi- 
tions may be significantly different in both frames. We will discuss 
this issue in more detail in a forthcoming article (Gozdziewski et 
al., in preparation). 

Each run of the hybrid code has been initiated by random se- 
lection of the GA population (between 512 and 4096 individuals), 
considering possibly wide parameter ranges. For instance, the range 
of orbital periods was set blindly to [800, 63600] days, and angles 
and eccentricities were set to their whole possible ranges. The orig- 
inal "population" was then transformed by GA operators over 512- 
1024 generations. Each member of the final set was then used as 
an initial condition for the L-M algorithm, and the resulting solu- 
tions were sorted and stored. The hybrid procedure was repeated 
hundreds of times for each combination of model-data set. We ex- 
amined whether the obtained solutions converged to the same min- 
ima. Due to the semi-deterministic nature of the GAs, one should 
interpret the results in a statistical sense. 

The same procedure may be applied to the Newtonian model, 
as the planetary contribution x can be computed independently of 
the optimisation method. (In this case the derivatives to the LM 
algorithm were approximated numerically). 

Finally, uncertainties of the best-fit parameters were deter- 
mined using the bootstrap algorithm (Press 2002), as the variances 
of parameters in a tested solution that has been re-fitted to 4096 syn- 
thetic data sets drawn randomly with replacement from the original 
sample. We found that due to the particular distribution of OPTIMA 
observations that are grouped in small "clumps" of a few data 
points, the bootstrap algorithm tends to underestimate the uncer- 
tainties when compared to the formal error determination through 
the diagonal elements of (Xv)^^^ curvature (covariance) matrix. 



3 KINEMATIC MODELING THE LITERATURE DATA 

To verify the literature models of the HU Aqr system, we gathered 
Barycentric Julian Dated (BJD) egress times published by Schwope 
et al. (2001), Schwarz et al. (2009) and Qian et al. (201 1). That data 
set consists of 1 13 points, and will be called the SSQ set hereafter. 
Our first attempt was to reproduce the results of Qian et al. (201 1) 
with our formulation of the LTT model. We did not expect this to 
be straightforward, since their model assumes the inner planet to be 
on a circular orbit. We conducted calculations for two ephemeris 
models, linear and quadratic (Eqs. 10 and 11), respectively. 



3.1 The linear kinematic ephemeris 2-planet model 

In the linear ephemeris case, we found many, almost equally good 
2-planet solutions with (Xy)^^^ ~ 1.15 and an rms ~ 2.3 sec. In 
these best fits the inner planet has a period of ~ 5500 — 6000 days. 
However, the period of the outer planet varies between 7000- 
20000 days. The resulting systems imply (0-C) residuals caused in- 
dividually by the planets in wide ranges, up to ^ 6000 seconds, and 
companions in basically any mass, eccentricity and period range 
while still preserving excellent rms ~ 2.3 sec and similar "flat" be- 
haviour of the residuals. The left panel of Fig. 2 shows the most 
exotic and actually the best-fit solution found in our experiment. 
The Keplerian fit parameters of this solution, as well as its inferred 
elements are given in Table 1 (Fit A). This fit is very different from 
those found by Qian et al. (2011), Wittenmyer et al. (2012), and 
even in the last paper by Hinse et al. (2012). This configuration 
has (Xv)^^^ ^ 1.143 and an rms ~ 2.3 sec, and is characterised by 
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almost equal orbital periods of ~ 5470 days. The pericenter argu- 
ments of the planets in this fit differ by nominal value of 180.2° 
and as a result, the Keplerian barycentric orbits are almost exactly 
anti-aligned, with planets placed close to their periastrons at the 
initial epoch. This configuration could be understood as a pair of 
Trojan-planets in Ic: lb mean motion resonance (MMR). Although 
the resulting LTT signal has apparently small amplitude ~ 60 sec- 
onds as shown in the (0-C) diagram (see the left-hand panel in 
Fig. 2), the LTT semi-amplitudes ^b,c are excessively large (up to 
~ 6000 seconds), implying just absurdly massive companions of 
^ 10, 000 Jupiter mass each (~10 M© !). This solution reveals that 
an inherent degeneracy of the LTT signal (and its model) may ap- 
pear because the signal is the result of the differential gravitational 
tugs of the companions on the binary. Indeed, in this particular Tro- 
jan configuration, even small deviations from the anti- alignment of 
orbits leads to large changes in the planetary masses (over 3 orders 
of magnitude) and semi-major axes (within a range of a few AU), 
indicating that as they are not supported by the currently available 
observations, these dynamical parameters are poorly constrained. 
The mathematical fit permits putative companions as massive as 
stars but in reality, such objects should influence dynamical and 
spectral properties of the binary system. Such solutions are there- 
fore excluded. 

The lc:lb MMR solution is a vivid example demonstrating 
that due to the possibility of configurations involved in extremely 
strong mutual interactions, modeling the LTT signal globally (with- 
out any a priori assumptions on the system configuration) cannot 
be studied in terms of the kinematic model. In general, an exact, 
self-consistent A/^-body model should be used to determine the ini- 
tial conditions. 



3.2 The quadratic kinematic ephemeris 2-planet model 

In the case of a quadratic ephemeris, we found a well defined min- 
imum of (Xv)^^^ ^ 0.972, which is an apparently statistically per- 
fect solution. Its synthetic curve with measurements over-plotted 
is shown in the right-hand panel of Fig. 2, and orbital parameters 
are given in Table 1 as Fit B. That solution has been frequently 
found in different runs of the hybrid code, which reinforces its 
global character. To show the latter, we computed parameter scans 
of (Xv)^'^^ the (Pb,c,^b,c)-planes (Fig. 3), by fixing points of a 
grid in a given plane and minimizing (Xv)^^^ ^^^^ all remaining 
free parameters of the model. This made it possible to obtain stan- 
dard confidence levels as marked with coloured curves. The best-fit 
solution is again very different from solutions found in the litera- 
ture. While the elements of the inner planet are well constrained, 
the orbit of the outermost companion reveals extremely large ec- 
centricity (~ 1). That points again to highly degenerate (unrealistic) 
best-fit solutions, with near-parabolic or even hyperbolic, open or- 
bit of one "planet" — as the fit implies — being a low-mass stellar 
object of ~ 160 Jupiter-masses. Other solutions with slightly worse 
(Xv)^^^^ 1-1 and still very similar rms ~ 2.2 sec may be found too, 
which means that the quadratic ephemeris model is unconstrained 
by the SSQ data. 

In the quadratic ephemeris model, the orbital periods are close 
to the 4:3 ratio, which is equivalent to the low-order 4c :3b mean 
motion resonance. In addition, the eccentricity of the outer planet is 
extreme, close to 1 . Hence again, the kinematic formulation seems 
inadequate to derive the proper initial condition of the multiple- 
planet configuration. We conclude here that when we only have the 
SSQ data at our disposal, there seems to be no unique and phys- 



Table 1. Jacobian geometric parameters with the inferred masses and semi- 
major axes of 2-planet LTT fit models in the barycenter frame with the 
linear and parabolic ephemeris to the SSQ data set. Synthetic curves with 
data sets are shown in two panels of Fig. 2 and (Xv)^^^ scans in Fig. 3. 
Numbers in parentheses represent the uncertainties at the last significant 
digit. Total mass of the binary is 1.08 Mq (Schwarz et al. 2009). Indices 
"b" and "c" refer here to planets "1" and "2" in the mathematical model 
Eqs. 10 and 1 1. See the text for more details. 



Model 


Fit A 


FitB 


parameter 


linear ephemeris 


parabolic ephemeris 






10.2 ± 0.5 


Pi frlnv*;! 

r]j LCiaysj 


JH-U / HZ H- i 


9Q1 f) -1- 98 




0.138 ±0.034 


0.34 ± 0.07 


co^ [degrees] 


207 ± 21 


22 ±8 


n [BJD 2,440,000+] 


11694 ± 175 


5652 ± 97 


Kc [seconds] 


5942 ± 17 


322 ± 30 


Pc [days] 


5476 ± 424 


3931 ±50 


ec 


0.141 ±0.035 


0.99 ± 0.05 


cOc [degrees] 


27 ±20 


358.3 ± 0.2 


Tc [BJD 2,440,000+] 


6214 ±451 


9207 ± 42 


Pbin [days] 


0.086820400(4) 


0.0868204250(8) 


To [BJD 2,440,000+] 


9102.92004(2) 


9102.91988(2) 


(3 [X 10-13 cycle-2] 




-3.06(6) 


at [au] 


1.375 


4.08 


Mb sin / [Mjup] 


9780 


5.69 


ac [au] 


1.374 


4.58 


mcsini [Mjup] 


9811 


159 




113 


113 




1.143 


0.972 


rms [seconds] 


2.31 


2.13 
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Figure 4. OPTIMA hexagonal fibre bundle centred on HU Aqr. The ring 
fibres (1-6) are used to monitor the background sky simultaneously. 



ically meaningful solution explaining the LTT variability. Or, the 
planetary fit model and its assumptions are incorrect. 



4 NEW OBSERVATIONS AND DATA REDUCTION 

4.1 Observations with OPTIMA and other instruments 

To resolve the model degeneracies as described above, we gathered 
new, yet unpublished observations of the HU Aqr binary. The new 
collected mid-egress BJD times are given in Table 2. These data ex- 
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HU Aqr 2-planet LTT-fit (linear ephemeris, SSQ data), Vxr=1 -143, rms=2.3 s 



HU Aqr 2-planet LTT-fit (quadratic ephemeris, SSQ data), Vxr=0.972, rms=2.13 s 




= 5467 ±41 8 days 6^ = 0.14 ±0.04 
: 5476 ± 424 days 6^ = 0.14 ±0.04 
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Figure 2. Synthetic curves of the best-fit, 2-planet solutions to the mid-egress BJD times of the SSQ data set. The left panel corresponds to a Hnear ephemeris 
model (Eq. 10) and the right panel corresponds to a quadratic (parabolic) ephemeris model (Eq. 11, the parabolic trend has been removed). Panels are labeled 
with the orbital periods and eccentricities of the putative companions. Bottom plots with shaded background show the residuals after subtracting planetary and 
astrophysical contributions from the LTT signal. Discontinuous-like features of the parabolic ephemeris model around cycles and 46,000 appear due to the 
extreme eccentricity of the outer body. See Table 1 for the orbital and inferred elements (Fit A and B, respectively). 
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Figure 3. Parameter scans of the (Xv)^^^ function computed around the best fit solution to the SSQ data for the quadratic ephemeris model (see the right-hand 
panel in Fig. 2 and Fit B in Table 1). Colour curves are for the formal l,2,3a-levels of the best-fit solution whose parameters are marked with an asterisk. 



tend the work of Schwope et al. (2001), Schwarz et al. (2009) and 
Qian et al. (201 1). The currently available data set of HU Aqr egress 
times consists of 171 measurements in total, including 10 points 
presented in Qian et al. (2011). Among these measurements, 68 
were obtained with the Optical Timing Analyzer (OPTIMA) in- 
strument that operates mostly at the 1.3 m telescope at Skinakas 
Observatory, Crete, Greece. 

The high-speed photometer OPTIMA is a sensitive, portable 
detector to observe extremely faint optical pulsars and other highly 
variable astrophysical sources. The detector contains eight fibre- 
fed single photon counters — avalanche photodiodes (APDs), and 
a GPS for the time control. There are seven fibres in the bundle 
(Fig. 4) and one separate fibre located at a distance of ~ l^ Sin- 
gle photons are recorded simultaneously and separately in all chan- 
nels with absolute UTC time-scale tagging accuracy of ~ 4 ^s. 
The quantum efficiency of the APDs reaches a maximum of 60% 
at 750 nm and lies above 20% in the range 450-950 nm (Kan- 
bach et al. 2003, 2008). During the HU Aqr observations, OP- 
TIMA was pointed at RA(J2000) = 2lh07"^58.n9, Dec(J2000) = 
— 05°17'40f'5, corresponding to the central aperture of the fibres 
bundle (Fig. 4). For sky background monitoring, we usually choose 
one out of the six hexagonally located fibres. We look for the fibre 



that is not by chance pointed to any source, therefore records sky 
background, and its response is the most similar to the central fibre 
response when the instrument is targeted at the dark sky. An exam- 
ple of a sky background subtracted light curve is shown at the top 
left-hand panel in Fig. 5. 

We derived new fits to the HU Aqr eclipse egress times, as 
well as reanalysed many of the already published OPTIMA data. 
There are 26 eclipses obtained by OPTIMA published by Schwarz 
et al. (2009). We were able to reanalyse only 21 out of the 26 
light curves, because only those were available in the OPTIMA 
archive. Our completely new data set includes 42 precision photo- 
metric observations, starting from cycle / ~ 29, 900, overlapping in 
time window with the literature data. We derived 23 new eclipse 
profiles from the OPTIMA data archive spanning 1999-2007 and 
obtained 19 new OPTIMA optical HU Aqr light curves in 2008- 
2010. Note that some of the OPTIMA observations have been al- 
ready published in the very recent literature (Nasiroglu et al. 2010). 

We also gathered and reduced 1 1 observations performed at 
the MONET (MOnitoring NEtwork of Telescopes) project which is 
a network of two 1.2 m telescopes operated by the Georg-August- 
Universitat Gottingen, the McDonald Observatory, and the South 
African Astronomical Observatory. These precision data in white 
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light (500-800 nm) were binned in 5 s intervals, with 10~^ days 
(0.1 s) accuracy, separated by 3 s readout. The most recent obser- 
vations were performed at November 18th, 201 1. 

An additional three egress times were obtained from the 
eclipse observations carried out with the PIRATE telescope 
equipped with the SBIG STLIOOIE CCD camera (Holmes et al. 
2011). PIRATE, funded by the Open University Department of 
Physics and Astronomy, is a remote-controlled telescope located 
at the Astronomical Observatory of Mallorca (OAM), Spain. 

We also performed optical observations of HU Aqr in white 
light with the 1.5-m Carlos Sanchez Telescope (TCS) equipped 
with Wide FastCam (WFC, Fig. 6). The WFC is a Ik x Ik - pixel 
camera with optics offering a FOV of 12 arcmin with a scale of 0.6 
arcsec/pix. HU Aqr eclipses were observed on September 30, and 
October 01, 201 1 with integration times of 3 and 5 seconds, respec- 
tively. WFC works in frame transfer mode, therefore the readout 
time is effectively null or in other words is equivalent to the expo- 
sure time. UTC mid-exposure times of the photometric measure- 
ments were converted to the Barycentric Julian Dates in Barycen- 
tric Dynamical Time using the procedure developed by Eastman 
et al. (2010). 

Some technical details of the observations performed with the 
MONET/N, TCS, and PIRATE telescopes are given in Table 3. 



4.2 Determining time markers of the eclipses 

In Fig. 5 we show an example of HU Aqr high time resolution 
OPTIMA photometric (see the right top panel, and blue curve in 
the left-hand top panel) and polarimetric (Stokes I, red curve in the 
right-hand top panel) light curves. These graphs are to be compared 
with the light curves from TCS, obtained with 3 and 5 seconds 
exposures illustrated in Fig. 6. Obviously, the OPTIMA resolution 
makes it possible to track the egress phase closely, which enabled 
us to determine the mid-egress moments very precisely. 

Measuring the time of mid-egress properly is critical to ob- 
tain the (O-C) diagrams, since it is the time marker of the eclipse 
(Schwope et al. 2001; Schwarz et al. 2009). To derive the mid- 
egress moment ^o, the sigmoid function 

parametrised by a 1,^2 and was fitted to the light curve points 
in the egress phase of the eclipse, spanning preselected exponential 
scaling parameters . We found that there is no strong dependence 
of the derived to on the adopted At. This can be seen in the bot- 
tom left-hand panel of Fig. 5 where three mid-egress times are 
marked with black open squares. These moments are derived for 
three different choices of A^: 0.1, 0.5 and 1.0, respectively. While 
these times depend on At, they fall within a 2 second range, as 
marked by a shaded strip at the bottom left-hand panel of Fig. 5. A 
half of that range (^ 1 second) may be typically estimated as the 
maximum possible error of to in the OPTIMA data set. The formal 
lo uncertainty of the sigmoid fit in this case is still smaller and at 
the level of ~ 0. 1 second. Moreover, the shape of the eclipse may 
significantly depend on the spectral window. Panels in the right col- 
umn of Fig. 5 illustrate the light curves of HU Aqr derived in the 
optical, white band domain (the blue curve) and in the polarimet- 
ric domain (Stokes I, the red curve). In the latter case, the egress 
looks quite different and spans over a longer time. Given that these 
two data sets were taken four years apart, the observed difference 
might have been caused by different emission states of the source. 



Table 2. New HU Aqr BJD mid-egress times on the basis of light curves 
collected with: OPT-ES022 — OPTIMA photometer installed at ESQ 
(Chile), OPT-SKO — OPTIMA operated at the Skinakas Observatory 
(Crete), PIRATE — a telescope at the Astronomical Obs. of Mallorca, 
MONET/N — the network of telescopes at the McDonald Obs. and the 
SAO (South Africa), and WFC — the 1.5-m TCS (Canary Islands). 



Cycle 


BJD 


Error [days] 


Instrument 


29946 




0000037 


OPT-FS022 


29957 


945170^ 7QQ^545 


000003 R 




29958 


945170^ RR61705 


0000034 


OPT-FS022 


30265 


24517^0 5400324 


0000041 


OPT-SKO 


30287 


2451732 4500Q02 


0000023 

\J ,\J\J\J\J\J J 


OPT-SKO 
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0000033 


OPT-SKO 


30300 
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0000054 
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2451734 ddfiQldO 


0000031 


OPT-SKO 
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2451734 5337R56 


000001 R 


OPT-SKO 
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0000030 
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00000R4 
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000001 5 

yj .yjyjyjyjyj i j 
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0000024 

yj .yjyjyjyjyj 


OPT- SAO 

y J I I vj r\y J 


44534 


2452969 3R00760 

Z^^^Z^yKjy .uOKJKJ 1 KJKJ 


0000033 

yj .yjyjyjyjyj J J 


OPT-NOT 


44557 


2452Q71 376Q377 


00000R5 

yj .yjyjyjyjyj J 


OPT-NOT 


51020 


2453532 4Q71 5Q5 


0000100 

\J .\J\J\J\J i \J\J 


OPT-SKO 
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2453536 4Q0Q030 


0000064 

yj . yjyjyjyjyjyj^ 


OPT-SKO 

WJT X vjXVW 


51067 


2453536 '^11121^ 


0000033 

yj . yjyjyjyjyj J j 


OPT-SKO 

WJT X OXVW 


55535 


2453Q24 4Q1 3426 


0000102 

yj . yjyjyjyj x wz^ 


OPT-SKO 

WJT X vJXvW 


55627 


2453Q32 47RR164 


0000061 

yj .yjyjyjyjyjyj i 


OPT-SKO 

WX X LJXVW 


55661 


2453Q35 4307071 


0000064 

yj .yjyjyjyjyjyj^ 


OPT-SKO 


55719 


2453940 4662754 


0000162 

yj .yjyjyjyj i vjz 


OPT-SKO 
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24543 IQ 524240Q 


0000074 
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OPT-SKO 
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jl^'-tu^ I L\j.'-T\j 1 L'-ryyj 
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\j.\j\j\j\j\jj J 
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OPT-SKO 
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0000032 
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OPT-SKO 
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OPT-SKO 
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MONET/N 
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MONET/N 
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MONET/N 


77031 


2455790.7824571 
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MONET/N 


77066 
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0.0000077 


MONET/N 


77067 


2455793.9079841 


0.0000055 


MONET/N 


77078 


2455794.8630179 


0.0000065 


MONET/N 


77546 


2455835.4949490 


0.0000179 


WFC 


77557 


2455836.4499905 


0.0000295 


WFC 


77789 


2455856.5922852 


0.0000038 


MONET/N 


77802 


2455857.7209399 


0.0000090 


MONET/N 


77823 
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0.0000066 
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78100 
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Figure 5. High time resolved OPTIIVIA light curves of HU Aqr. The upper left panel: photometric eclipse in the white light. The bottom left panel: a close-up 
around the mid-egress overplotted with three fitted sigmoid functions with different values set for the parameter (0.1, 0.5, 1.0, respectively). The fitted mid- 
egress times are denoted by open squares. The shaded region corresponds to a time span of 2 seconds. The upper right panel: a comparison of photometric and 
polarimetric OPTIIVIA light curves of HU Aqr obtained in 2004 and 2008, respectively. The bottom right panel: a close-up of photometric and polarimetric 
HU Aqr Hght curves around the mid-egress time showing the difference between the egress shapes. The shaded region covers 6 seconds. 



Table 3. Technical data of the IVIONET/N, PIRATE and TCS observations 
of HU Aqr. Tobs represents the observation time span, and AT is the mean 
exposure time of a single frame. 



Date Instrument Filter Tobs [hours] AT [s] 



2010 Oct 06 


PIRATE 


WL 


1 


10 


2010 Oct 19 


PIRATE 


WL 


3.5 


10, 20 


2011 Sep 30 


TCSAVFC 


WL 


1 


3 


2011 Oct 01 


TCSAVFC 


WL 


1 


5 


2011 Apr 03 


IMONET/N 


WL 


0.60 


8 


2011 May 03 


IVIONET/N 


WL 


0.50 


8 


2011 Jul 21 


IMONET/N 


WL 


0.32 


8 


2011 Aug 17 


IVIONET/N 


WL 


0.60 


8 


2011 Aug 20 


IVIONET/N 


WL 


0.46 


8 


2011 Aug 20 


IVIONET/N 


WL 


0.10 


8 


2011 Aug 21 


IVIONET/N 


WL 


0.58 


8 


2011 Oct 22 


IVIONET/N 


WL 


0.25 


8 


2011 Oct 23 


IVIONET/N 


WL 


0.50 


8 


2011 Oct 25 


IVIONET/N 


WL 


0.16 


8 


2011 Nov 18 


IVIONET/N 


WL 


0.45 


8 



To derive the mid-egress moments gathered in Table 2, the sigmoid 
function was fitted to all single light curves. 



4.3 On the light curves in different spectral windows 

Vogel et al. (2008) and Vogel (2008) obtained high time resolved 
and accurate light curves of HU Aqr during its low state using the 
ULTRACAM (Beard et al. 2002; Dhillon et al. 2007) at the Very 
Large Telescope (VLT, May 13, 2005). These authors decomposed 
the light curve into three emission components emerging from the 
accretion spot, the photosphere surrounding it, and from the white 
dwarf itself. As a result, they were able to derive the temperature of 
the WD ~ 13500(200) K and the temperature of the accretion spot 
~ 25500(1500) K. They also estimated the ratio of the spot area to 
the WD surface to be on the level of 5%. The black body spectra of 
the WD and of the spot have their maxima at 215 nm and 113 nm, 
respectively. The accretion spot and the accretion stream are time 
variable in brightness, as well as in the geometric position in the 
system. Therefore, the orbital phase at which they occur is not con- 
stant. The ULTRACAM delivers simultaneous light curves in three 
colours: w, g and r. An example of such a three colour HU Aqr light 
curve can be found in Fig. 6 of Schwarz et al. (2009), where the 
shape — energy dependence can be easily seen. Thus, a comparison 
of egress times in different wavelength domains was possible. In 
the two cases with filters u and r, the WD constitutes the main con- 
tribution to the egress intensity, because it can be seen unperturbed 
when it comes out of eclipse. However, during high and intermedi- 
ate accretion states, the WD might be out-shined by the accretion 
stream. In the w-band, the spot contributes 25% of the emission. 
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Figure 6. HU Aqr eclipse light curves obtained with the WFC mounted on the TCS. Observations were performed on 30/09/2011 and 01/10/2011 with 
integration times of 3 and 5 seconds, respectively. Bottom panels are for the close-up of eclipse egress, with fitted sigmoid function (solid line, see Eq. 12). 
Blue vertical lines mark the determined mid-egress times. 



while in the r-band its contribution is only 12% (Vogel et al. 2008; 
Vogel 2008). This suggests that as the time marker of the eclipse, 
it is better to use more "reddish" than "blueish" data, particularly 
in our case as we have broadband observations gathered in X-rays, 
UV and optical domains at our disposal. 

There exists evidence that the EUVE light curves differ from 
quasi- simultaneous ROSAT/HRI light curves, as it can be seen, for 
example, for eclipses recorded on October 1996 and May 1997, as 
shown in Fig. 2 of Schwope et al. (2001). The eclipse ingress is 
often not measured because of strong suppression of soft X-rays 
by absorbing matter along the accretion stream. When the eclipse 
duration can be determined, the eclipse duration seems shorter in 
the case of EUVE data. Thus the derived mid-egress moments can 
be shifted by a few seconds. According to Schwope et al. (2001), 
an expected variation of the eclipse span should be not more than 
0.001 of the orbital phase which corresponds to not more than ~ 8 
seconds. Also Schwope et al. (2004) show in their Fig. 3 evidence 
of a different eclipse length as well as phase folded egress shapes 
at soft X-rays, HST UV, and high-speed optical photometry with 
a multichannel multicolor photometer (the MCCP 2.2 m telescope 
at Calar Alto) during the 1993 high state and the 1996 low accre- 
tion state. The scatter of the egress times resulting from changes 
of the accretion geometry during high and intermediate accretion 
states is estimated on the level of 2 seconds (Schwarz et al. 2009). 
Large differences between light-curves due to the eclipse of the 
accretion stream are also visible in the optical photometric mea- 
surements performed in parallel with the ROSAT observations (see 
Fig. 3 in Schwope et al. 2001). Some of those light curves were 
obtained with rather poor time resolution, e.g. 53 and 12 seconds. 

It is worth mentioning that time stamps calculated by Schwope 
et al. (2001) and Schwarz et al. (2009) for the photon counting UV 
and X-ray detectors were computed from the mean of the arrival 



times of the first three photons after the eclipse, while for the opti- 
cal observations, they used the moments of the egress half intensity, 
which is common in the literature. Examples of the HU Aqr light 
curves obtained by the XMM-Newton EPIC-PN and Optical Mon- 
itor detectors are presented in Figures 2, 3 and 4 of Schwarz et al. 
(2009). XMM observations, contrary to the bright state ROSAT ob- 
servations, were not resolved at time scales shorter than 2 seconds 
due to the low count rates. 

We first used all available egress times, archival as well as 
new ones, to model the (0-C) diagram. However, given the above- 
mentioned arguments, we decided to select only those measure- 
ments that were obtained in the white light or photometric V band, 
in order to keep the data more uniform and homogeneous. This ap- 
proach renders the measurements independent of possible varying 
emission regions in different bands. We also decided to skip the 
most "suspicious" egress-times at some stage of fitting the orbital 
model, which is described below. 

We note that the HST observations (three points around / ^ 
14,000) were performed with the EOS instrument in the 120- 
250 nm range. These points were also excluded in our further anal- 
ysis, falling out of the white light and the V band range. 



5 MODELING ALL RECENT DATA 

Thanks to the new set of precision OPTIMA mid-egress measure- 
ments, as well as observations performed at PIRATE, TCS and 
MONET/N telescopes, we can re-fit planetary models to the whole 
set of data up to November 18th, 2011. We fitted the data with the 
linear and quadratic ephemeris models (Eqs. 10, 11). 



© 0000 RAS, MNRAS 000, 000-000 



On the HU Aquarii planetary system hypothesis 1 1 



HU Aqr 1 -planet LTT-fit (linear ephemerls, all data), ^^=^QA, rms=12.6 s 



35 
25 

¥ 15 

c 

o 

o c; 

0) 

9 -5 

o 

-15 
-25 





1 

Pb = 227^ 


1 

5±20d 


ays 
i 


1 1 1 
L T eb = 0.99±0 

^Wi ' 


1 1 

01 




1 1 


1 1 1 


* — 

i 

1 1 * 


1 1 1 1 1 1 1 


1 1 1 1 1 1 1 



10 20 30 40 50 60 70 

cycle number ['lO^] 



Figure 7. Synthetic curve of the 1 -planet LTT model with linear ephemeris 
to all available data, including the very recent egress times collected by the 
OPTIMA photometer, as well as PIRATE, TCS and MONET/N telescopes. 
Open circles are for measurements in Qian et al. (201 1). 



5.1 Single-planet models to all recent data 

At the first attempt, we tested the 1 -planet hypothesis. For the linear 
ephemeris model, the 1 -planet solution is characterized by extreme 
eccentricity and displays large residuals and a strong trend present 
in the (0-C) diagram (see Fig. 7). This suggests a more general 
quadratic model, on which we focus now. 

The results derived for the whole set of 171 measurements are 
shown in the top panels of Fig. 8. Interestingly, the 1 -planet model 
fits the data very well in a large part of the time- window between 
/ = 25 , 000 and / = 80, 000 (see the left-hand panel of Fig. 8). How- 
ever, over approximately one fourth of the time-window (/ = to 
/ = 25, 000), the data fit the synthetic curve poorly. That can be bet- 
ter seen in the close-up of the residuals shown in the top right panel 
of Fig. 8. It appears that the residuals follow a regular and char- 
acteristic "damping" trend, that could be associated with a mass- 
transfer process ongoing in the binary or solar-like magnetic cy- 
cles. Results of our experiments show that the recent observations 
by Qian et al. (201 1) appear to be outliers to our 1 -planet solution, 
as the mid-egress times are shifted by about of 3-10 seconds w.r.t. 
the synthetic curve. Because these observations overlap in the time 
window with much more precise OPTIMA data, that discrepancy 
between these two data sets cannot be avoided. Actually, obser- 
vations by Qian et al. (2011) do not fit any model that has been 
tested with the OPTIMA observations, including 2-planet models 
and both types of the ephemeris (see the Appendix). 

In an effort to explain the strange behaviour of the residuals, 
we realized, as it was discussed already, that the available observa- 
tions come from different telescopes/instrumentation, and to make 
the matter worse, the egress times are measured on the basis of 
light curves in different spectral windows. In particular, the first 
part of the data set contains the egress times derived from X-rays 
(ROSAT and XMM) and ultraviolet (EUVE, XMM 0M-UVM2 
and HST/FOS) light curves, and some eclipses were observed with 
OPTIMA in polarimetric mode. To remove the possible inconsis- 
tency due to the different spectral windows and filters, we consid- 
ered data sets consisting of the egress times measured only in the 
optical range (white light and the V band). The results are shown 
in the bottom panels of Fig. 8 for the optical data without X-ray 
and UV, but including polarimetric measurements (note, that the 
polarimetry was done in the white-light band; compare with the 



Table 4. Keplerian parameters for the 1 -planet LTT fit model with quadratic 
ephemeris to all data gathered in this work (Fit I) and to measurements se- 
lected in the optical and V-band domain (Fit II). Synthetic curves with mid- 
egress times overplotted are shown in Figs. 8 and 9. Numbers in parentheses 
represent the uncertainty at the last significant digit. Total mass of the binary 
is 0.98 Mq (Schwope et al. 201 1). See the text for more details. 



Model 


Fit I 


Fit II 


parameter 


all measurements 


optical measurements 


[seconds] 


13.9 ±0.3 


14.7 ± 0.2 


Ph [days] 


3278 ± 28 


3287 ± 19 


eh 


0.03 ± 0.04 


0.13 ±0.04 


cOb [degrees] 


211 ±40 


226 ± 10 


7b [BID 2,440,000-h] 


6233 ± 360 


6361 ± 102 


Pbin [days] 


0.0868204226(5) 


0.0868204259(4) 


7b [BID 2,440,000-h] 


9102.92004(2) 


9102.91994(1) 


(3 [X 10-13 day cycle'^] 


-2.61 (5) 


-2.95(4) 


ah [au] 


4.29 


4.30 


mbsin/ [Mjup] 


6.71 


7.10 


A^data 


171 


115 




5.23 


2.48 


rms [seconds] 


4.8 


3.7 



top panels of Fig. 8 for all data gathered). As can be seen from the 
bottom panels in Fig. 8, the "damping" effect has almost vanished, 
suggesting that it could have appeared due to the presence of X-ray 
and UV-derived eclipses. Still, there is a group of data points with 
large errors, around / ~ 14,000, which do not fit well to the clear 
quasi- sinusoidal variation of the (0-C). The deviations of these 
points may be explained by poor time -resolution (~ 12 seconds of 
the AIP07 CCD camera), that has been used to observe the HU Aqr 
eclipses (Schwope et al. 2001). Let us also note that the Qian et al. 
(2011) data points are again systematically outliers with respect to 
the synthetic signal. After removing these data and all points (seven 
measurements) in the polarimetric mode, we obtained a homoge- 
neous optical data set to which we fitted the quadratic ephemeris 1- 
planet model again. The synthetic curve of this fit with data points 
over-plotted is shown in Fig. 9. Parameters of this fit are presented 
in Table 4 as the final solution Fit II and are well constrained by 
the observations. To demonstrate the latter, we show projections of 
(Xv)^^^ in selected two-dimensional parameter-planes (see Fig. 10) 
close to the best-fit model. As can be seen, there is a strong corre- 
lation between the time and argument of pericentre which can be 
understood noting that the orbital phase (kh = G5b + Mb ) must be 
preserved. 

The best-fit model seem to constrain the damping factor P ~ 
—3 X 10~^^ day cycle"^ very well. Such a value is close to es- 
timates in the literature, e.g., ^ — 5 x 10~^^ day cycle"^ by 
Schwarz et al. (2009) and - -2.5 x 10"^^ day cycle"^ by Qian 
et al. (2011). It is still larger by more than one order of mag- 
nitude to be explained by gravitational radiation, but remains in 
the range of magnetic braking (Schwarz et al. 2009). A similar 
large-magnitude period decrease has been found in other CVs, like 
NN Ser (- -6 x 10"^^ day cycle"^, Brinkworth et al. 2006). Be- 
sides the angular momentum loss, the large magnitude of the period 
change is commonly explained as due to the Applegate mechanism 
(basically excluded in the HU Aqr) and/or a the presence of a very 
distant, long-period companion body. Likely, a few astrophysical 
and/or dynamical effects may be involved that could determine ap- 
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Figure 8. The top row: synthetic curve of the 1 -planet LTT model with quadratic ephemeris to all available data, gathered in this work, including the very 
recent mid-egress times collected by the OPTIMA photometer, as well as PIRATE, TCS and MONET/N telescopes (the top left) with orbital parameters given 
in Table 4 (Fit I), and close-up of residuals to that model {the top right panel). The bottom row: the same for the white light and visual band (V) data, including 
polarimetric observations by OPTIMA (i.e., the UV- and X-band observations are excluded) shown at the bottom left panel, and its residuals (the bottom right 
panel). The white filled circles mark the Qian et al. (201 1) measurements. 



HU Aqr 1 -planet LTT-fit (quadratic eph., optical w/o polarim.), A/xr=2.59, rms=4.1 s 
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Figure 9. Synthetic curve of the 1 -planet LTT model with quadratic ephemeris to observations in white light -i- V band (see the text for more details) without 
measurements in Qian et. al (2011) and polarimetric data. Orbital parameters of this solution are given in the Table 4 (Fit II). 
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Its synthetic curve with data points overplotted is shown in Fig. 9. The large symbol marks nominal elements of the solution. Closed curves are for formal 
1,2,3a levels of (x?)^^^ with scale displayed at the colour legend. 



parently secular period decrease. Its definite explanation is com- 
plex, and we consider this as a subject of a new, forthcoming work. 

We also fitted the quadratic ephemeris model only to the high- 
est precision OPTIMA data. The results for measurements that in- 
clude polarimetric observations are shown in Fig. 11. For that case, 
we found a period similar to the quadratic ephemeris model for the 
entire data set. The fit has very small rms ^0.8 sec. The relatively 
large (Xv)^^^ ~ 3.4 of the OPTIMA solution in this case may sug- 
gest that the adopted uncertainties, at the ~ 0.1-0.2 second level 
(in a large sub- set of the measurements) are in fact underestimated. 
We also identified the most deviating points as coming from the 
polarimetric measurements (see e.g., a point marked in the resid- 
uals plot around / = 65,000, and the residuals of both solutions). 
To examine, whether these data may change the solution, we fitted 
the quadratic ephemeris model to the white light OPTIMA obser- 
vations only, skipping all polarimetric data. The best-fit orbital pe- 
riod of ~ 3400 days remains close to the full-coverage window fit. 
A slightly smaller rms of ~ 0.7 seconds suggest a better fit with- 
out the polarimetric data, indeed. The orbital periods coincidence 
cannot be fully proved due to the relatively narrow observational 
window of the OPTIMA white light measurements. Actually, the 
parameter scan (not shown here) reveals that the la contour around 
the (x?)^^^ minimum in the ( Pb, ^b)-plane is "opened" on the right 
side of the orbital period-axis, hence it can not be constrained yet 
by the OPTIMA data alone. 

5.2 Alternate models to all recent data 

Finally, using the hybrid optimiser, we performed additional ex- 
periments by fitting three models to all available data: the 1 -planet 
model with a heuristic sine damping term, and 2-planet models, 
both in terms of the linear and quadratic ephemeris. We also per- 



formed N -body modeling of the 2-planet configurations. The re- 
sults, which are described in the Appendix, imply that all these 
models lead to non-unique solutions or configurations with simi- 
lar orbital periods, for which the kinematic model is inadequate, as 
we discussed above. Some of these kinematic best-fit 2-planet so- 
lutions are qualitatively similar to configurations found for the SSQ 
data set, with orbital periods ratios close to Ic: lb and 4c :3b, respec- 
tively. The extended data set still does not constrain the Keplerian 
2-planet models. 

The same can be concluded for the A/^-body models (see Ap- 
pendix, Sect. A3). Although we found stable configurations in 
terms of the quadratic ephemeris, the semi-major axis of the outer 
companion is unconstrained (between 4 AU and at least 20 AU). 
These stable fits exhibit varying sign of P (it means that the binary 
period might decrease or increase). Moreover, stable solutions with 
relativelty small p may be found in very narrow stability zones in 
the (^c^c) -plane, see the Appendix and Fig. A6. These areas are 
associated with low order MMRs, like 3c :2b MMR. It is very uncer- 
tain though, how massive companions of HU Aqr could be locked 
in such tiny stability areas. Hence, some larger values of provid- 
ing extended zones of stable motions seems more likely (see panels 
of stable fits labeled by IV, V, and VI of Fig. A6). However, there is 
also a correlation between the magnitude of P and the semi-major 
axis of the outer planet. For relatively distant planet c, |P| may be 
~ 2 X 10~^^ day cycle"^, which is difficult to explain by magnetic 
braking or mass-loss. However, this may indicate a presence of a 
third companion in an unconstrained orbit. 

We conclude that these results seem to favour the 1 -planet 
hypothesis as the simplest model explaining the (0-C) variability, 
particularly in the light of very small rms of the homogeneous OP- 
TIMA set and apparently perfect quasi- sinusoidal fit illustrated in 
Fig. 11. 



© 0000 RAS, IMNRAS 000, 000-000 



14 K. Gozdziewski et ah 



HU Aqr 1 -planet LTT-fit (quadratic ephemeris, all OPTIMA), Vxr=3.37, rms=0.79 s 



-15 





1 1 1 
.^-^ Pb = 3205 ± 396 days 
\ 6^ = 0.06 ±0.04 


1 1 


1 1 1 1 1 




1 1 1 

1 1 1 


1 s 1 

polarimetric point 




\ I \ i 

1 1 



20 30 40 50 60 70 

cycle number ['lO^] 



Figure 11. Synthetic curves of the 1 -planet LTT quadratic ephemeris mod- 
els to optical OPTIMA measurements, including polarimetric data One of 
the most deviating polarimetric points is labeled in the residuals panel. 



6 RED NOISE AND/OR SYSTEMATIC ERRORS? 

Analysis of the LTT observations has much in common with pul- 
sar timing, planetary transits, and precision radial velocity observa- 
tions, which are modelled with least-squares under the assumption 
that the measurement errors are uncorrelated (white noise). How- 
ever, as is known particularly by pulsar observers, the assumption 
that white noise is the only source of error is unjustified when aim- 
ing at estimating the underlying model parameters and their uncer- 
tainties (Coles et al. 201 1). In the past, this effect had been respon- 
sible for false detections of planets around pulsars (Bailes et al. 
1991). Similarly, correlated (red) noise or systematic errors have 
been found in the planetary transit data (Pont et al. 2006) and very 
recently, in the radial velocity measurements (Baluev 2011). The 
same type of non-Gaussian, low-frequency correlation of residuals 
to the orbital period of the binary may be present in the LTT data 
collected over long time intervals. 

The danger of such systematic effects in the LTT-analysed bi- 
naries is reinforced due to their activity and complex astrophysical 
phenomena responsible for the observed emission. One of the al- 
ready well recognized mechanisms able to produce cyclic variation 
of the orbital period of the binary has been proposed by Applegate 
(1992). As shown by this author, a magnetic star (here, the sec- 
ondary) changes its internal structure due to magnetic cycles. The 
latter implies a variable zonal harmonic coefficient J 2 and subse- 
quently, a variable gravitational tidal field for the orbital compan- 
ion which results in a varying orbital period (Hilditch 2001). The 
Applegate mechanism as a possible origin of large (0-C) variations 
in the HU Aqr data was studied in detail by Vogel (2008), as well 
as by Schwarz et al. (2009). They discarded this possibility since 
the HU Aqr stellar setup does not provide enough energy to drive 
changes of the orbital period. Similar results were obtained for the 
NN Ser system, that likewise has a low-mass, low-luminosity sec- 
ondary star (Brinkworth et al. 2006) with a conclusion that it is 
incapable of driving significant period changes in terms of the Ap- 
plegate model. 

Another mechanism explaining observed long-term periodic- 
ities could be a slow precession of the rapidly spinning magnetic 
WD star, which has been proposed as a source of long periods de- 
tected in a few CVs, for instance FS Aur and V455 And (Tovmas- 
sian et al. 2007). However, HU Aqr is unlikely to host such a WD, 
as this AM Her-like system is known to be synchronously locked. 



As a first, yet preliminary attempt, we tried to determine the char- 
acteristic that can be used to quantify the shapes of the HU Aqr 
light curves and might help to detect their variability and hence as- 
trophysical sources of the LTT residuals. This approach mimics the 
bisector velocity span (BVS) technique used to detect distortion of 
spectral lines due to stellar spots and chromospheric activity. It is 
well known that stellar spots may produce apparent radial velocity 
changes up to 200 ms ~^ (Berdyugina 2005). As a similar charac- 
teristic to the BVS, we choose the slope of the linear function fitted 
to the egress phase of the light curve, usually spanning no more than 
a few seconds interval. We analysed 59 available light curves in the 
precision OPTIMA set. The results are shown in Fig. 12. In seven 
cases, we decided the data were not precise enough to derive the 
slope reliably (as indicated by green filled squares) because of, for 
instance, bad weather or strong wind that could introduce telescope 
guidance errors. In a few other cases (seven again, marked with 
blue triangles), only two points were taken for the fit, and there- 
fore no error estimation was possible. Nevertheless, the obtained 
slopes are uniform and span less than 2 degrees range close to 90 
degrees. That furthermore indicates a similarly rapid egress phase. 
The results of this test are encouraging, and support the planetary 
hypothesis. 

However, the slopes should be best re-computed for all avail- 
able light curves that were used to determine the egress-times. The 
problem of the non-homogeneity of the collected light curves still 
exists. Due to varying eclipse profiles (e.g. during different accre- 
tion states), the determination of mid-egress dates is often very dif- 
ficult. For instance, it could be prone to rather subjective choices 
of the photometric data range to fit the parameters of the sigmoid 
function, Eq. 12. That may introduce significant systematic errors, 
particularly if the reduction is performed by different researchers. 
This issue may be likely resolved by a re-analysis of the entire set 
of all available light curves, under similar conditions paying partic- 
ular attention to their origin - the spectral window, an instrument, 
and even technical and observational circumstances. 

Another direction still open is a study of the binary interac- 
tions, to eventually eliminate or discover astrophysical causes of the 
LTT variability. The problem is in fact universal and affects other 
techniques of extrasolar planets detection, such as pulsar timing 
and radial velocity monitoring of active or evolved stars, as well. 
It is yet possible that the observed (0-C) signal has both the plan- 
etary and unmodeled astrophysical component (Potter et al. 201 1), 
making its unique resolution even harder. 

To the best of our knowledge, possible effects of the red-noise 
regarding the LTT observations have not been studied in detail. 
That problem certainly deserves a deep and careful investigation. 



7 CONCLUSIONS 

Using a new formulation of the LTT model of the (0-C) to the avail- 
able data of the HU Aqr system, we found that the 2-planet hypoth- 
esis by Qian et al. (2011) is not likely. Our results reinforce recent 
negative tests of dynamical stability of that system in the literature. 
The self-consistent LTT model presented in this work exhibits de- 
generate solutions, such as (apparently) Trojan objects of ~ 10"^ 
Jupiter masses, or a companion in an collisional/open parabolic 
or hyperbolic orbit. Ironically, two such solutions to the literature 
SSQ data are the best-fit models found in extensive, quasi-global 
searches adopting a hybrid optimization. 

Moreover, on the basis of a much extended, precision data set, 
collected by the OPTIMA network, that increased the number of 
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Figure 12. Linear slopes of 59 HU Aqr egresses derived on the basis of OPTIMA light curves by fitting a linear function only to the egress phase, usually 
spanning not more than a few seconds. In seven cases, the light curves were not precise enough (as indicated by green filled squares) because of a bad weather. 
In other seven cases, only two points were taken for the fit, therefore no error estimation was possible (blue triangles). See the text for more details. 



data points analysed in previous works by ~ 50 %, we have shown 
that the observed (0-C) variations may be consistently explained 
by the presence of only one circumbinary planet of the minimal 
mass of ~ 7 Jupiter-masses, in an orbit with a small eccentricity 
of ^0.1 and an orbital period of ~ 10 years, similar to Jupiter 
in the Solar system. Our results support the original 1 -planet hy- 
pothesis by Schwarz et al. (2009) rather than the 2-planet model 
proposed by Qian et al. (2011). If confirmed, that planet would be 
the next circumbinary object detected from the ground, shortly af- 
ter such companions have been announced around HW Vir (Lee 
et al. 2009), NN Ser (Beuermann et al. 2010), UZ For (Potter et al. 
2011), SZ Her (Lee et al. 2011), DP Leo (Beuermann et al. 2011), 
followed by recent discoveries of Kepler- 16b (Doyle 201 1), Kepler- 
34b and Kepler-35b planets (Welsh et al. 2012). According to esti- 
mates by these authors, the observed rate of circumbinary planets 
around close binaries may be ~ 1%. 

Also, we found that the observations by Qian et al. (201 1) are 
not confirmed by the OPTIMA measurements due to systematic 
relative shift of ~ 3-10 seconds. The nature of this discrepancy is 
yet unknown. If the shift is caused by an error, all 2-planet models 
presented in the literature that make use of their data are affected. 

Besides the disagreement between our conclusions and the 
previous works, our results suggest that the kinematic modeling of 
2-planet configurations is not fully justified on the grounds of the 
dynamics because the best-fit models may imply large masses (up 
to stellar range), large eccentricities, and similar orbital periods in- 
dicating a possibility of strong mean motion resonances. Moreover, 
the (0-C) variability that suggests 2-planet solutions most likely 
appears due to mixing observations done in different spectral win- 
dows. That feature of the data set - as we have shown here - intro- 
duces systematic effects that may alter the best-fit solutions signif- 
icantly. This conclusion is supported by extensive numerical sim- 
ulations of the 2-planet systems dynamics by Horner et al. (2011), 
Wittenmyer et al. (2012) and Hinse et al. (2012). Considered within 
statistical error ranges, the initial conditions lead to catastrophically 
disruptive configurations, unconstrained elements of the outermost 
body, and/or period damping factor |3. 

In this work, we found best-fit stable 2-planet models within 
the quadratic ephemeris N-bo&y model to all available data, but 
the semi-major axis of the outer planet cannot be yet constrained. 
Stable configurations are located within low-order MMRs spanning 
tiny stable zones in the phase space, or are characterised by a large 
magnitude of the period decrease. In the first case, it is difficult 
to explain, how a few Jupiter mass companions could be trapped 



in such particular, isolated resonances. In the second case, a large 
I PI requires an efficient, internal mechanism of the binary period 
change, or indicates a presence of one more companion. Our find- 
ings might be a breakthrough after a few cited works reporting ba- 
sically only unstable 2-planet models of the HU Aqr system, but 
these discrepancies add even more ambiguity to the 2-planet hy- 
pothesis. 

However, the results of our experiments show that the 1 -planet 
solution is relatively well constrained by available optical observa- 
tions selected as a homogeneous data set. Because the early opti- 
cal data (the white light and V-band measurements) are coherent 
with an impressive, very clear quasi- sinusoidal signal exhibited by 
superior-precision OPTIMA measurements, as well as with the re- 
cent MONET/N, PIRATE and WFC data, a single-companion hy- 
pothesis seems well justified. A confirmation of the planetary origin 
of the LTT signal still requires long-term monitoring of the system. 
Due to its very long orbital period, it will take many years to con- 
firm or reject the signal coherence. Such new data would be also 
very useful to constrain the orbital period by the recent OPTIMA 
observations alone. 
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APPENDIX A: ALTERNATE MODELS TO ALL DATA 

In this section, we display supplementary results illustrating a few 
alternative models to the 1 -planet solution of the (0-C) of the HU 
Aqr binary that was analysed in the main part of the paper. Ba- 
sically, all 171 data points are modeled, although in some cases, 
we removed outlying data from Qian et al. (2011) because they 
clearly introduce a systematic error. We considered 1 -planet model 
with a heuristic, sine-damping term (Sect. Al), 2-planet kinematic 
models (Sect. A2) and the full, 2-planet, self-consistent Newto- 
nian model (Sect. A3). The aim of this Appendix is to demon- 
strate that 2-planet models lead to non-unique or unconstrained 
solutions. Hence, these results reinforce a hypothesis of a single, 
quasi- sinusoidal signal of possibly planetary origin. 



Al Quadratic ephemeris 1-planet model with damping term 

To describe the suggested damping signal visible in Fig. 8 (the top 
right-hand panel), we modified the quadratic ephemeris model by 
adding a heuristic term having the following form: 

TdampCO =''^0+Aexp(-^/rdamp)sin(ndamp^ + C|)o), (Al) 

where Tq is an offset, A is the semi-ampHtude of the signal, Tdamp is 
the damping time scale, ^damp = 27c/Pdamp is the frequency, and (|)o is 
the initial phase at epoch / = 0. Two examples of best-fit solutions 
to all available data (171 measurements) are shown in Fig. Al. Let 
us note that the planetary orbital period in the configuration in the 
right-hand panel of this figure is twice the period in the 1-planet 
model studied earlier. The (0-C) of a solution shown in the left- 
hand panel cannot be distinguished from 2-planet models (see the 
text below). 

A physical nature of that damped signal is uncertain. Allowing 
for some speculations, the damping might appear due to a long- 
term relaxation in the binary system which may, for instance, be 
due to the binary's magnetic cycles (Applegate mechanism). In 
such a case, the observed LTT signal would be resulted from two 
distinct phenomena. However, we recall here that Vogel (2008) and 
Schwarz et al. (2009) estimated that the Applegate mechanism can- 
not be responsible for orbital changes of HU Aqr. 



A2 Kinematic 2-planet models 

We also tested 2-planet models with the linear and quadratic 
ephemeris, (Eqs. 10, 11). Examples of the best-fit configurations 
with comparable (Xv)^^^ ™s are shown in Fig. A2. Similar 

to the case of the SSQ data set, no unique solution may be found. 
For the paraboHc ephemeris, we found many similar-quality best- 
fit solutions. These fits are characterised by the orbital periods ratio 
close to lc:lb MMR with inferred planetary masses of ~ 20 M 
jup (the bottom left-hand panel in Fig. A2), close to 4c: 3b MMR 
with inferred masses of ~ 5.5 M jup and ~ 4.0 Mjup (the top right- 
hand panel in Fig. A2). A solution close to the 2c: lb MMR (the 
top right-hand panel in Fig. A2), as well as configurations with ex- 
treme eccentricity Cc ^ 0.95, positive damping factor P ^ 10~^^ cy- 
cles day~^, and unconstrained Pc ^ 250,000 days (not shown here) 
was also found. In all these cases, the rms remains at the level of 
2.4 seconds. Some of these solutions are qualitatively similar to the 
2-planet fits found for the SSQ data set. These results imply that 
the significantly extended data set still does not constrain 2-planet 
models. 

For the linear ephemeris model, we found one best-fit solu- 
tion that frequently appeared in different runs of the hybrid code. 
It is shown in the bottom right-hand panel of Fig. A2. This solu- 
tion is characterised by an orbital periods ratio close to 6:5. Taken 
literally, this fit corresponds to Trojan brown dwarfs. However, the 
kinematic model is inadequate for such a configuration of massive 
objects. 



A3 Newtonian, self-consistent A/^-body 2-planet models 

In the light of the discussion presented above, we performed a pre- 
liminary modelling of all available data with the help of the hybrid 
algorithm driven by the self-consistent A^-body model. Moreover, 
we tested Lagrange stability of the best-fit models following their 
orbital evolution over at least 10^ orbital periods of the outermost 
planet. Configurations which survived during such time without a 
collision or remaining on closed orbits were regarded stable. In this 
experiment we use 161 data points, excluding data in (Qian et al. 
2011), due to the discrepancy with OPTIMA measurements. 

To illustrate the results of the hybrid optimization, we pro- 
jected the found solutions onto particular planes of the Keplerian 
astrocentric, osculating elements of the planets (Fig. A3) at the 
epoch of the first observation. The general finding is that the N- 
body formulation helps to improve the rms, that decreased from 
~ 2.4 seconds to ~ 1.9 seconds as compared to kinematic models. 

The top row of Fig. A3 illustrates the results for the linear 
ephemeris. Clearly, the data do not constrain the semi-major axes 
and eccentricities of the companions. The eccentricities tend to be 
large, up to 0.8. Moreover, the best-fit configuration exhibit simi- 
lar values of semi-major axes (^5.6 AU and ^6.3 AU) and large 
masses in the brown-dwarfs range of ~ 20 Jupiter masses. We did 
not find any stable configurations within this model. It is consistent 
with the results for the SSQ data set (Hinse et al. 2012). 

Interesting results are obtained for the quadratic ephemeris 
model (see the bottom row in Fig. A2) although also this model 
does not constrain orbital parameters, due to even larger spread of 
the semi-major axes and eccentricities than in the linear ephemeris 
model. Two minima of (Xv)^^^ ^^e found, around a\3 ^ 4 au, and 
flb ~ 6 au, respectively. The best-fit configurations have (Xv)^^^ 
2.6 and and an rms ~ 1.9 second that is ^ 20% better than for 
the best kinematic models. In the neighborhood of the first (Xv)^^^ 
minimum ( a^^ 4 AU), we found a few thousands of Lagrange 
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HU Aqr 1 -planet LTT-fit (quad, ephem. + damped sine, all data), V%^=2.6, rms=2.4 s 



HU Aqr 1 -planet LTT-fit (quad, ephem. + damped sine, all data), Vx^=3.6, rms=2.6 s 
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Figure Al. Examples of synthetic curves of the 1 -planet LTT quadratic ephemeris models with a sine damping term to all 171 data points gathered in this 
work, including 10 measurements in Qian et al. (201 1), marked with white filled circles. The shaded curves represent the planetary and the damped sine signal, 
respectively. 



HU Aqr 2-planet LTT-fit (quadratic ephemeris, all data), V%r=2.41 , rms=2.35 s 



HU Aqr 2-planet LTT-fit (quadratic ephemeris, all data), Vx^=2.42, rms=2.35 s 
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Figure A2. Examples of synthetic curves of the 2-planet LTT quadratic and linear {the bottom right panel) ephemeris models to all 171 data points analysed in 
this work, including 10 data points in Qian et al. (201 1) which are marked with white filled circles. The shaded curves represent single planetary signal terms, 
respectively. 



stable models characterised by (Xv)^^^ < ^ and an rms < 2.1 (still 
better than for the best 2-planet kinematic models). These fits have 
well bounded ab ~ 4 AU and small eccentricities up to 0.4. How- 
ever, the osculating semi-major axis of the outer body is uncon- 
strained and covers many low order mean motion resonances, be- 
tween 3c:2b MMR and 5c: lb MMR. Figure A4 shows synthetic 
curves of two example solutions corresponding to the 3c:2b MMR 



(the left panel) and for a model close to 3c: lb MMR (the right 
panel). To identify these resonances, in the neighborhoods of a 
few selected best fit models, we derived high-resolution dynamical 
maps (1440 x900 data points) shown in Fig. A5. These maps are 
computed in terms of the fast indicator MEGNO (Cincotta et al. 
2003), with the help of our recently developed CPU cluster soft- 
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Figure A6. Statistics of 2-planet A/^-body quadratic ephemeris models gath- 
ered with the hybrid algorithm, projected onto the (<3c , P)-plane, for the 
quadratic ephemeris. Meaning of symbols is the same as in Fig. A3. 

ware MECHANIC^ (Slonina et al. 2012). MEGNO measures the 
maximal Lyapunov exponent, which makes it possible to distin- 
guish between chaotic and regular solutions. Each point at these 
maps has been integrated over ~ 10^ orbital periods of the outer- 
most companion. The dynamical maps confirm that the Lagrange 
stable models examined over a limited time- span are equivalent to 
quasi-periodic, stable solutions. 

A solution illustrated in the left panel of Fig. A4 is the best- 
fit stable model found in the hybrid search with (Xv)^^^ ^ 2. 82 
and an rms~ 2 sec. It is located in a very narrow, isolated stabil- 
ity island of the 3c:2b MMR and characterised by relatively large 
P ^ —2.7 X 10~^^ day cycle"^, similarly to the kinematic model. 
The right panel of Fig. A4 shows a configuration close to the 
3c: lb MMR, which has even larger p ~ — 6 x 10~^^ day cycle"^. 

Figure A6 shows a statistics of the best-fit solutions in the 
(^c, P)-plane. It reveals that P is not constrained, regarding even 
its sign. Stable models exhibit a strong correlation between both 
these parameters. A larger value of the semi-major axis of the outer 
planet is related to a larger magnitude of p. Due to this correlation, 
an interpretation of stable configurations is complex. For relatively 
small magnitude of p, stable configurations are characterized by 
low order MMRs and may be found in tiny areas of stable motions 
(see Fig. A5). For more separated planets, when stability zones 
are much more extended, |P| increases. Already |P| ~ 2 x 10~^^ 
day cycle"^ is difficult to explain by physical phenomena in the 
binary, as we discussed in Sect. 5. Such large values of |p| may in- 
dicate a third, long-period companion object in a very distant orbit. 
However, because already the 2-planet model is not constrained by 
the data, also a 3 -planet configuration cannot be fixed without am- 
biguity. We did an attempt to search for such Newtonian 3 -planet 
models within the linear ephemeris, but we did not find any im- 
proved, nor stable solutions of this type. 
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Figure A3. Statistics of 2-planet A/^-body models gathered with the hybrid algorithm, projected onto planes of selected parameters. The top row illustrates 
the results for the linear ephemeris, the bottom row shows the fits for the quadratic ephemeris model. The rms quality of these solutions is coded as filled 
circles: the better fit — the darker colour. Orbital parameters of the best- fit solutions are marked with shaded, intersecting lines and a flower symbol. Solutions 
Lagrange-stable over 10^ revolutions of the outermost planet are marked with red- white circles. 
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Figure A4. Examples of synthetic LTT curves and (O-C) residuals for stable 2-planet quadratic ephemeris Newtonian (A^-body) models to 161 data points 
analysed in this work, without 10 data points in Qian et al. (2011). These models are selected from a sample illustrated in Fig. A3. Dynamical maps of these 
solutions shows Fig. A5 (they are labelled with the Roman numerals, as I and IV, respectively). 
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Figure A5. MEGNO dynamical maps in the {uc , ec)-plane for a few representative A/^-body stable solutions illustrated in the bottom panels of Fig. A4. Yellow 
colour encodes strongly unstable (chaotic) configurations, and purple colour (MEGNO (F) ~ 2) is for stable, quasi-periodic solutions. Parameters of the 
nominal, tested fits are marked with the star symbol. The most prominent, low-order mean motion resonances are labeled. The original resolution of these 
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